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Abstract 


This research investigates residual-based a posteriori error estimates for finite element 
approximations of heat conduction in single-layer and multi-layered materials. The 
finite element approximation is based upon hierarchical modelling combined with 
p-version finite elements. Hierarchical modelling results in thermal elements which 
are geometrically compatible with structural finite elements. Thermal stresses are an 
important concern when designing reusable launch vehicles, and accurate temperature 
distributions throughout the structure are required. 

A posteriori error estimation is a way to determine the accuracy of an approxi- 
mate thermal solution when the exact solution is unknown. Error estimates are also 
necessary for developing an adaptive scheme which is the automatic process of mesh 
refinement and p-enrichment to deliver a solution with the desired accuracy. 

Element error indicators are determined by solving an element equation for the 
error using the element residual, and a global error estimate in the energy norm is 
computed by collecting the element contributions. Two methods, the average flux and 
the equilibrated flux method, are discussed for constructing the element flux boundary 
condition for the error equation. The error estimation is extended to multi-layered 
materials, and a directional error indicator is developed to distinguish the error in the 
hierarchical model from the error in the finite element method. Comparisons of the 
actual and estimated error show that the equilibrated flux method provides accurate 
estimates of the error for single and multi-layered materials. Numerical results also 
show that the directional indicators accurately determine which contribution to the 
total error dominates. This is an essential step in implementing an effective adaptive 


scheme. 
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Chapter 1 
Introduction 


When performing structural analysis of hot structures, it is necessary to understand 
the temperature distribution throughout the structure. Thermal concerns include 
maximum and minimum temperatures, thermal stresses induced by temperature gra- 
dients, and rates of heat flow depending on the application of the material being stud- 
ied. Therefore, efficient and dependable methods of computing thermal solutions are 
vital in the design and analysis of thermal structures. For example, reusable launch 
vehicles must be able to repeatedly withstand the extreme temperatures present dur- 
ing Earth re-entry. The Thermal Protection System (TPS) is the primary outer 
material that protects the system from high temperatures. Materials and concepts 
are continuously being investigated to improve the insulating properties for TPS. 
Such multi-layered materials include composite laminates and TPS panels composed 
of insulation between two layers of metal to replace the tiles used on the orbiter of the 
Space Shuttle system. Therefore, the numerical analysis of these types of materials, 
along with error estimates, will aid in the design process. 

1.1 Background 

The finite element method is a popular numerical analysis tool for engineers working 
thermal problems with complex geometry or boundary conditions. The advances in 


1 



finite element approaches and computer technology have increased the usefulness of 
finite element methods and allowed for more accurate and efficient approximations. 
One such advanced method is the hierarchical p- version finite element method. The 
use of hierarchical modelling allows the dimension of the domain to be reduced by 
one dimension. The benefit of this approach is that mesh regeneration is not required 
when changing the approximation order to obtain a more accurate solution. 

When making numerical approximations with the finite element method, it is 
useful to know how accurately the methods approximate the exact solution. Error 
estimation is widely known and used to answer the question, ’How good is the finite 
element solution?’ Knowledge of the application of the problem and an understanding 
of how the error behaves can reduce the computational effort needed to achieve the 
desired accuracy. 

The two types of error estimates are a priori and a posteriori . A priori error 
estimates contain the exact solution and the parameters that influence the accuracy 
of the approximate solution. This type of error estimate provides information on 
the convergence and stability of the method and gives the asymptotic behavior of 
the error in the approximation as mesh size and polynomial orders are appropriately 
varied [1]. 

The second type is a posteriori error estimates, that are computed from the finite 
element approximation and the given problem data and do not require the exact 
solution to the problem. Since approximation methods are generally used because the 
exact solution is not known for a the problem of interest, a posteriori error estimates 
are useful for determining accuracy when analyzing complex problems for which the 
solution is unknown. A posteriori error estimates are also necessary for adaptivity 
and control of the finite element approximation error. Adaptivity is a method of 
automatically refining the mesh or increasing the polynomial degree in particular 
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regions where the finite element error is greatest. These regions of large error occur 
near steep solution gradients, singularities, or discontinuities in applied loadings. The 
mesh can be refined or the polynomial degree increased in these particular areas 
to improve the solution instead of the costly approach of refining the entire mesh. 
Through the use of a tolerance, local refinements can be made after each approximate 
solution is obtained until the error is within a specified tolerance, and this repeating 
method can be done with no further input from the analyst. Therefore the method 
can be set up and allowed to continue without any further input from the analyst. 
Accurate a posteriori error estimates are essential to the success of such an adaptive 
procedure. 

The two main categories of a posteriori error estimates are explicit and implicit. 
Explicit error estimates involve a direct computation (usually post-processing) using 
the finite element solution and the given problem data and include an unknown con- 
stant which is typically ignored. Implicit error estimators involve the approximation 
of a boundary value problem for the error using residuals and generally involve the 
solution of a system of equations, and are usually more accurate than explicit error 
estimates [2]. Although more computational time is generally required for implicit 
error estimates, the improved accuracy can be worth the extra effort. 

When a mesh is defined for finite element analysis, the process usually involves a 
combination of experience, intuition, and guesswork. If the results of the finite element 
approximation appear reasonable, then the solution is accepted. If not, then the mesh 
is redesigned. But this process is time consuming, and without an a posteriori error 
estimate, there is no other reliable and precise way of judging the acceptability of the 
solution except to uniformly refine the mesh until the solution converges. A posteriori 
error estimation provides the analyst with a method of measuring the quality of the 
computed solution. A more accurate solution can then be obtained by selecting 
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various parameters to change such as the value of the polynomial approximation 
and the mesh size locally where errors are large. For the hierarchical p- version finite 
element method, directional error indicators indicate which parameters will efficiently 
improve the solution by providing a measure of the modelling error and the finite 
element error separately. 

1.2 Review of Previous Research 

The hierarchical p- version finite element method was pioneered by Babuska, Szabo, 
and Peano in the mid to late 1970’s and early 1980’s [3] and [4]. Hierarchical mod- 
elling helps to simplify the problem by reducing the dimension by one before the 
p-version finite element method is applied. An optimal set of basis functions for 
the through-thickness direction that account for the piecewise continuous behavior 
of the solution to multi-layered problems was developed by Vogelius and Babuska in 
1981 [5]. Another method for approximating the solution in multi-layered problems 
is the Zig-Zag method developed in 1996 by Averill and Yip [6], which uses a single 
polynomial approximation through the thickness with a piecewise linear function su- 
perimposed on it. The coefficients for the piecewise linear function are determined 
from continuity conditions at the interfaces of the layers. The optimal basis functions 
are used with the hierarchical p-version finite element method for multi-layered mate- 
rials in this research since the use of the Zig-Zag method requires more computation 
when using p-enrichment for the through thickness polynomial approximation. 

The subject of a posteriori error estimation has grown in popularity since the 
pioneering work of Babuska and Rheinholdt in the late 1970’s. The first a posteriori 
error estimates for linear elliptic problems were developed by Babuska and Rheinholdt 
to guide local mesh refinement [7] and obtain accurate results without refining the 
entire mesh. Since then, a posteriori error analysis has been developed for parabolic 
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and hyperbolic partial differential equations. One method of a posteriori error esti- 
mation is the flux-projection method, which is based on the smoothing of the fluxes 
and comparing them with the finite element fluxes [2]. The method used in this work, 
called the element residual method, was first introduced in 1984 at a conference in 
Lisbon [8] and uses the residuals of the finite element solution. The element resid- 
ual method requires element boundary conditions. One method of approximating 
the element boundary fluxes is the equilibrated flux method, which is based on the 
work by Ladeveze and Leguillon [9], Kelly [10], and Bank and Weiser [11] and is used 
in this research. Ainsworth and Oden advanced the element residual method using 
equilibrated fluxes for application to different types of problems [12]. 

The basic techniques for a posteriori error estimation were established in the early 
1990’s. In recent years, a new approach in the study of a posteriori error estimates has 
emerged called goal-oriented error estimation. The error is measured with respect to 
a specific goal of the analysis, called the quantity of interest, instead of in the energy 
norm, and techniques have been established to obtain upper and lower bounds of 
the error. This approach requires the solution of the adjoint problem, or the dual 
solution, to compute an influence function which relates the residual to the error 
in the quantity of interest [13]. Therefore, it is more computationally expensive. 
Currently, the emphasis is on the study of robustness of existing estimators and the 
identification of limits on their performance. 

1.3 Purpose 

The purpose of this research is to apply a posteriori error analysis to two-dimensional 
heat conduction for single and multi-layered materials. The three primary goals of 
this research are to (i) compare the results of a posteriori error estimates using two 
methods of approximating the interior boundary flux, an average flux and an equi- 
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librated flux, for single-layered materials, (ii) extend and examine the performance 
of the error estimate for multi-layered materials, and (iii) develop directional error 
indicators to distinguish between the modelling error and the finite element error. 

1.4 Scope 

The hierarchical p- version finite element method for solving steady state heat con- 
duction problems is presented in Chapter 2. The weak form of the boundary value 
problem is developed followed by the hierarchical modelling in the through-thickness 
direction and the spatial p-version finite element method. Basis functions are de- 
scribed, and element matrices are defined. Finally, the method for enforcing the 
boundary conditions and solving the global system of equations are discussed. 

The element residual method of a posteriori error estimation is discussed in Chap- 
ter 3. The derivation of the element boundary value problem for the error is presented. 
The element boundary conditions are discussed with two methods of approximating 
the interior element flux, the average flux and the equilibrated flux method. The 
solution for the error in the energy norm is presented next with some numerical 
examples. Finally directional error indicators are introduced with some numerical 
examples to verify the error estimates. Chapter 4 presents the application of element 
residual a posteriori error estimation to multi-layered materials. The finite element 
method for multi-layered materials is discussed with the choice of optimal through- 
thickness basis functions for hierarchical modelling. The boundary value problem for 
each layer is discussed followed by a discussion of the equilibrated flux method used 
to approximate the interior element boundary flux. Finally, some numerical examples 
are presented for a two layer problem. A review of the work presented in the thesis 
is discussed in Chapter 5. A summary of the first four chapters is presented followed 
by conclusions and recommendations for future work. 
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Chapter 2 

Finite Element Thermal Analysis 


In all but some simple cases, the exact solution of structural problems is complex 
and requires an approximate solution using numerical techniques. The finite element 
method is a powerful tool for numerically approximating the solution to a wide range 
of engineering problems. There are numerous methods available for constructing a 
finite element approximation. Traditional finite elements use linear approximations 
of field variables over each element. By increasing the number of elements using mesh 
refinement, the approximate solution converges to the exact solution. An alternative 
to this approach, the p- version finite element method, is to use higher-order elements, 
which assume a polynomial of order p for the approximation on each element. By 
increasing the value of p, or p-enrichment, the approximation can be improved without 
increasing the number of elements. The drawback to this method using conventional 
higher-order elements is that the mesh must be regenerated for increased values of p 
since the number of nodes in an element depends on p. Mesh refinement is also an 
option for the p- version method. 

A more practical method which has been developed is the hierarchical modelling 
combined with p-version finite elements method. This method uses elements with 
a fixed number of nodes by adding more basis functions to the existing nodes when 
the approximation order is increased. The use of hierarchical modelling allows the 
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reduction by one dimension of the domain which reduces the complexity. The benefit 
of this approach is that when the approximation order is changed, mesh regeneration 
is not required. This chapter describes the hierarchical modelling combined with p- 
version finite elements approach used for this research to approximate the solution 
for a heat conduction problem with specified heat flux and temperature boundary 
conditions. 

2.1 Boundary Value Problem for Conduction Heat 
Transfer 


The goal is to solve for the temperature, u, in a built-up structure using the general 
heat conduction equation. According to the conservation of energy, the rate of heat 
entering and leaving a body, the rate of heat generation inside the body, and the rate 
of heat stored by the body must be balanced 

PZv-Q t + V <7 = <? (2- 1 ) 

where q is the heat flux vector defined as [q x , q y , q z ] T , V is the gradient operator 
defined for a rectangular coordinate system as [J~, J“] , and p and c p are the 

density and the specific heat of the material, respectively [14]. The internal heat 
generation, Q, is a scalar function of position and is associated with conversion of 
some other form of energy, such as chemical, electrical, or nuclear, to thermal energy. 
The heat flux is related to the temperature gradient by Fourier’s Law 

q = —kVu (2.2) 


where n is the thermal conductivity tensor, which is defined for an anisotropic mate- 
rial as the following matrix of thermal conductivities 


k 

k 

k 


XX 

yx 

zx 


k k 

• v xy h^x 

k k 

rVyy TVy 

k Z y k z 


K = 


(2.3) 



Since heat flows in the direction of decreasing temperature, and the gradient points 
in the opposite direction, the minus sign is included in (2.2) to make the heat flow 
a positive quantity. Substitution of (2.2) into (2.1) results in the partial differential 
equation for temperature 

pcp^ - V T (kVu) = Q (2.4) 

To allow for comparison with exact solutions, the steady state heat conduction in 
a simple two-dimensional domain made of an orthotropic material will be considered 
in this research. The simplification to two dimensions from three dimensions is made 
by assuming the domain is infinitely long in the ^/-direction. The domain, Cl = 
(0, L) x (jy, is shown in Figure 2.1, and has length L and thickness d with the 
boundary denoted by dCl. The thermal conductivity of a two-dimensional orthotropic 



material is different in the two principle directions, in this case the x and z directions, 
and is assumed to be constant. 


* 0 


0 k z 


(2.5) 


For steady-state conditions, (2.4) becomes a second-order elliptic partial differential 
equation. For this research, the boundary conditions are defined on each edge as either 
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a Dirichlet or a Neumann boundary condition. A Dirichlet boundary, denoted by T D , 
is defined as a boundary on which the temperature is specified and is also called an 
essential boundary condition. A Neumann boundary, denoted by IV, is defined as a 
boundary on which the heat flux is specified and is also called a natural boundary 
condition. The outward unit normal to the boundary is denoted by uq = [n x , 0, n z ] T . 
The boundary value problem for the temperature, u(x, z ), in the domain is given by 


V T (kVu) 

Q 

in 


(kVu) t tiq 

= Qs 

on Tjv 


u 


on T d 

(2.6) 


where q s and u s are not necessarily constant. 

2.2 Weak Formulation 


The finite element method is based on the weak form, or the variational form, of the 
boundary value problem (2.6). A weak form is a weighted- integral statement of a 
differential equation in which the differentiation is distributed among the dependent 
variable and a weight function and includes the natural boundary conditions of the 
problem [15]. To obtain the weak form of the boundary value problem for heat 
conduction, (2.6) is multiplied by a scalar test function v and integrated by parts 


[-v T (nVu)v<m = 

[ Qv dQ, 

(2.7) 

Jn 

Jn 

1 + [ (kVu) t Vv dQ = 

\dn J Q v ' 

[ Qv dfi 
Jn 

(2.8) 

+ [ (kVjj) t Vv = 

Tjv Jn 

[ Qv dfi 
Jn 

(2.9) 


The test function is required to be zero on T D , which simplifies (2.9) to a general 
problem statement for the temperature in the domain where the boundary conditions 
in (2.6) have been applied. Find the temperature u = u(x, z) such that u = u s on 
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and 


[ (nVufVvdQ = [ Qvdto- [ q s v ds (2.10) 

Jq Jq JT} v 

for all admissible test functions v where s denotes a boundary. If B(u, v ) and L(v) 
are defined as 


B(u,v) 


(kVu) t Vv dti 


m 



( 2 . 11 ) 


then the problem statement can be expressed in abstract form. Find u E y(fl) such 
that u = u s on T^ and 


B{u , v) = L{v) for all admissible v E V(Q) (2.12) 

where V (fi) is defined as a subspace of H 1 (12), the Hilbert space consisting of functions 
with square integrable first derivatives [16]. 

2.3 Finite Element Method 


The solution to (2.12) is also the solution to (2.6) and is not an approximation. But 
since the solution space V(Q) is infinite-dimensional, it is difficult to solve for the 
exact solution. The finite element method is used to construct a finite-dimensional 
subspace of F(fi), denoted by V^fi), by subdividing the domain into a collection of 
elements and selecting a set of basis functions. The abstract form of the finite element 
method is to find u E V (J2) such that u = u s on T D and 

B(u,v) = L(v) (2.13) 

for all admissible test functions v E F(Q). This allows an approximation to the exact 
solution to be found, and the larger the choice for the subspace V^fi), the closer to 
f/(fl) it becomes, therefore making the approximation more accurate. The following 
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sections describe the method used to approximate the solution to (2.12), hierarchical 
modelling combined with p- version finite elements. 

2.3.1 Hierarchical Modelling 

In this research, hierarchical modelling refers to the process of describing the solu- 
tion on a domain whose spatial dimension is reduced by one. The geometry of the 
rectangular domain in Figure 2.1 can be described as a one-dimensional domain, a 
straight line of length L with a constant thickness d. More general geometries can 
be described in this manner by allowing the thickness to vary along the length. This 
approach is commonly used in structural mechanics to represent beams, flat plates, 
and curved shells. 

The approach used in structural mechanics to represent the two-dimensional so- 
lution on the one-dimensional domain is usually a global approach. The solution is 
assumed to have a polynomial distribution through the thickness of the domain, and 
the degree of the polynomial, denoted by p z , is assumed to be constant through- 
out the domain. This assumption transforms a single partial differential equation 
on the original domain into a system containing p z + 1 differential equations on the 
dimensionally-reduced domain. The finite element method is then applied to each 
differential equation in the system. In the hierarchical modelling approach used in 
this research, the solution is assumed to have a polynomial distribution through the 
thickness of the domain after the dimensionally-reduced domain has been subdivided 
into elements. This allows the hierarchical model order, p z , to vary along the length 
of an element, and therefore throughout the domain. Also, the hierarchical model 
order, p z , can be easily increased as part of the finite element method. 

For simplicity, the hierarchical modelling approach will be described for a two- 
node linear finite element. The extension to higher-order (p-version) finite elements 
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will be described in the next section. The one-dimensional domain representing the 
rectangular plate is divided into N elements as shown in Figure 2.2. As with tradi- 



tional finite element methods for one-dimensional problems, the solution on a typical 
element, fi e , is assumed to vary linearly along the length of the element, and is written 
in terms of a local coordinate system 

= 2 “ & u ° + 2 ^ Ul (2-14) 


The mapping for each element in the global domain from (x 0 , Xi) to the local domain 
(— 1, 1) is given by 


2{x — x 0 ) 

X\ — Xq 


( 2 . 15 ) 


However for hierarchical modelling, the nodal unknowns w 0 and U\ are not constants, 
they are functions of the through-thickness direction 

= q (1 — 0 u o( z ) + «(1 + £)ui( z ) ( 2 - 16 ) 
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where the subscripts 0 and 1 denote the left and right nodes of the element respec- 
tively. The unknown functions uq{z) and U\(z) are assumed to be polynomials of 
degree p z0 and p zl 

PzO 

Uo = i>i(v)uoi 

i= 0 

Ui =Y^tpi{v)uu (2.17) 

i=0 

where p zQ and p z i are positive integers, p = is a normalized through-thickness 
coordinate, ^ are polynomial basis functions, and Uji are the unknown constants, also 
called the element degrees of freedom. The complete expression for the approximate 
solution on an element is obtained by substituting (2.17) into (2.16) to get 

1 PzO 1 Pzl 

«ln* = «(! -05Z^i(*?)«oi+ (2-18) 

Z i = 0 Z i= 0 

The basis functions for the through-thickness polynomial approximations are the 
integrated Legendre polynomials, which are shown to be well suited for computer 
implementation and have favorable properties for numerical stability by Babuska et. 
al in [4] and are given by 

^o = l 
ih = V 

1 9 ? 1 pT] 

= y— -R-i (£)<*? 

= ~r===(Pi(v) ~ Pi- 2(r?)), * = 2,3,...,p, + l (2.19) 

V 2 (2*-l) 

where Pi is the Legendre polynomial of the first kind of order i given by [17] 

P 0 (v) = 1 
Pi(v) = V 

PM = 1 [(2i - l)vPi-M - (* - l)Pi- 2 iv)\ i = 2,3,... ( 2 . 20 ) 

L 
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The integration in (2.19) is accomplished by applying the properties of Legendre 
polynomials. Note that the usual linear approximation for a one-dimensional element 
is a special case of the hierarchical model, obtained when p z0 = p z ± = 0 in (2.18). 
When the degree of the through-thickness polynomial does not vary along the element 
length, that is when p z o = p z i = p z , (2.18) can be written as 


p * fi 1 

^loe = 9(1 - O U 0 i + ~(1 + Quu 

i = 0 L Z Z 

1 Pz 

j = 0 «=o 


fpiiv) 


( 2 . 21 ) 


2.3.2 P-Element Approximation 


The traditional approach to higher-order one-dimensional finite elements requires 
that an element with a polynomial approximation of degree p x contain p x + 1 evenly- 
spaced nodes. The basis functions for such an element are the Lagrange polynomials, 
and the element is referred to as a Lagrange element. A Lagrange polynomial basis 
function is associated with each node of the element along with a corresponding degree 
of freedom which is the value of the unknown at that point. Quadratic and cubic 
Lagrange elements and the corresponding basis functions are shown in Figure 2.3. 

The p- version approach to higher-order, one-dimensional, finite elements is to add 
a single node in the middle of the element, and a set of higher-order polynomials 
that are zero at the endpoints are added to the usual linear basis functions associ- 
ated with the end nodes. The number of basis functions, and hence the number of 
unknown coefficients, associated with the middle node is determined by the degree 
of the polynomial approximation p x . Specifically, there are p x — 1 basis functions 
associated with the middle node. Elements of this type are often referred to as hier- 
archical p-elements, because increasing the order of the elements is accomplished by 
simply adding a basis function to the existing set, as indicated in Figure 2.4, without 
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$=-1 ^- V , $= V , q =1 



(a) (b) 


Figure 2.3: One-dimensional (a) quadratic Lagrange element and (b) cubic Lagrange 
element with corresponding basis functions. 


the requirement of mesh regeneration. The finite element approximation in a typical 
element is written 

^ln e = 9(1 ~0 w o + ~(1 + 0 M i + = (2.22) 

Z Z i=2 i = 0 

The degrees of freedom corresponding to the middle node do not represent the value 
of the solution at the midpoint and are therefore denoted by the symbol a*. The 
unknowns u 0 and Ui are combined into the set of unknowns, a*. In other words, 
a 0 = Uo and a\ = U\. The basis functions used for the one-dimensional element, 
i.e., for the polynomial approximation in the x-direction, are similar to the through- 
thickness basis functions and are given by 

fa = ^(i-0 
<k = ^(i + 0 
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Figure 2.4: Element basis functions, <f> and for p = 1 through 7. 

<t>i = \ ^r- r 

= , 1 , i = 2,3, . . . ,p x + 1 (2.23) 

V 2 ( 2 * - 1) 

The procedure for combining hierarchical modelling with p- version elements is the 
same as described in the previous section. The constant unknowns in (2.22) are 
replaced by polynomial functions in the through-thickness direction 

"1 "1 p X 

«ln- = o (! “ 0«o(*) + « (! + 0«i(«) + £ ( 2 - 24 ) 

Z Z i=2 i=0 

where 


PzO 

U 0 = ^2ipi(r])uoi 
i = 0 
P*i 

«i = ' 52 ' t Pi(v)uu 

i = 0 

Pz2 

*j = (2-25) 

*=o 
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Again, if the hierarchical model order is constant for the element, p z0 = p z i = p z2 = p z , 
then the approximate solution in an element can be written 

Pz ^ ^ Px 

= E + + + fain) 

i = 0 [ 3 = 2 

Px Pz 

^ln e = ^ ( t ) i(C)' i Pi(v) 0j .n (2.26) 

j = 0 i=0 

The finite element approximation for the temperature distribution in an element is 
determined by the unknown coefficients, a^. The basis functions, and ^( 77 ) are 
polynomials of degree and p*, respectively. The number of unknown coefficients 
for the element is determined by (p x + 1 ) x (p z + 1 ). 

2.3.3 Element Matrices 


The weak form of the boundary value problem, (2.13), is an integral equation over the 
entire domain, which can be represented by a sum of integrals over each element in the 
domain. The domain of an element, denoted by fi e , is defined as (xq,X\) x (^, |), 
where Xo and X\ are the x-coordinates of the end nodes. The local contributions 
from the elements are combined into a global problem to allow for the solution to 
the finite element approximation. Note that an element end node is shared by a 
neighboring element. Since the unknown coefficients are associated with a node, the 
approximate temperature solution is continuous across element interfaces. Equation 
(2.13) is written as a sum of element contributions 


f (kVu) t Vv dn e = V f Qv d,Q, e - [ q s v ds 
Jn e , J n e Jdn e nr N 


provided that 

N e l 

= ft ( 2 . 28 ) 

e=l 

N e l 

< 90 e n Tjy = r N ( 2 . 29 ) 

e=l 
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Here N d denotes the total number of elements in the mesh and dQ. e D T N denotes 
the part of the boundary of an element which lies on IV. To simplify the derivation 
of the finite element matrices the approximate solution in (2.26) can be written in 
compact notation using outer tensor products 


Px Pz 

u\ne EE = (4>® = X T ae (2.30) 

i=0 , 7=0 

The symbol ® represents the outer tensor product, and <f> is a vector of length (p x + 1), 
^ is a vector of length (p z + 1), and % and a e are vectors of length (p x + 1) x (p z + 1). 
The outer tensor product [17] of an n x n matrix A with an m x m matrix B is given 
by 


i-l CN . . 

1 

A\2 * 
A.22 * 

* Ain 

* A2n 


' B n ■ 


— 

’ A n B 

A2lB 

A 12 B 

A 22 B 

• A\ n B 

■ A2nB 

A n \ 

A n 2 * 

A 


Bml 

Bmm 


A n \B 

A n 2B ■ 

■ • A nn B _ 
(2.E 


There are a variety of choices for the test function, ft, which results in a system of 
equations from (2.27). The Bubnov- Galerkin method, commonly referred to simply as 
the Galerkin method, defines the test function to be the same as the functions used in 
the approximation of the solution, ft, and will be used in this paper [18]. Substituting 
the approximate solution (2.30) for ft and XiA = 1, • • • , (p x + 1) x (Pz + 1) for ft, into 
(2.27) results in the system of equations 


Nel r -1 N e i 

E / Vx(«Vxf ffl e a e =5] 


e=l L' 


/ Qx dQ, e - [ q s x ds 

JYlr 


(2.32) 


Here the gradient of the vector x is defined as the gradient operating on each element 
of the vector. This equation can be written in matrix form as 


N e l N el 

E(A Q ea e ) = 5]TV 


e=l 


e=l 


(2.33) 
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The element matrix A ri e is defined in terms of a local mass and stiffness matrix using 
the outer tensor product as 


Afte — 




Vx(/tVxf dW 




dxdx T , dxdx T 


dx dx 


dz dz 


diY 


= k 


L 


X1 d<f) d<t> T i it j 

dx 0 V , V ,J dz + 


xq dx dx 


k z / dx 

Jxq 


■/ 


d ty 1 
f -| <9z dz 


dz 


k x K x ® + k z M x 0 JG 


(2.34) 


where M denotes a local mass matrix and if denotes a local stiffness matrix. The 
element load vector, E^e, is defined as 

F Q e = j Qx dCl e — [ q s x ds (2.35) 

Jn e J r e N 

2.3.4 Global System of Equations 

The global system equations, (2.33), are written 

Aa = F (2.36) 


where A and F are the global matrices formed by assembling each A^e and E^e, 
respectively, and a is the global vector of unknowns. The element vector of unknowns, 
a e , in (2.30), (2.32), and (2.33) is an element ordering of a with entries only for the 
unknowns associated with the nodes of the elements. The process of assembly is 
rather straightforward and can be found in any finite element textbook, such as [15], 
[19], and [20]. The approximate solution, u , is determined by applying boundary 
conditions to the global system of equations and solving the matrix equation for the 
unknowns a. 
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The Dirichlet boundary conditions in (2.10) are enforced on the system of equa- 
tions using the penalty method [19]. Recall the form of u from (2.26) 

Pz 1 1 Px 

= S o(! — O u oi + ^(1 + £)%* + ipiiv) (2.37) 

*= 0 |_ Z 3 = 2 

If the specified temperature is on the left side of an element, at £ = — 1, then <f)j (£) is 
zero for j = 2, ... ,p x and u becomes 

Pz 

«ln. (€ = -1) = T'MV)** (2-38) 

i= 0 

The value of u at this boundary is specified to be u s , therefore 

Pz 

(1)«00 + n fpi{v)uoi = u s (2.39) 

i= 1 

To enforce this boundary condition, we set Uoo = u s and Uoi = 0 for i = 1, . . . If 
the specified temperature is on the right side of an element, at £ = 1, we set U\q = u s 
and Uu = 0 for i = 1, . . . ,p*. Using the penalty method, the local coefficient matrix 
A^e and the local load vector F^e are modified by using a large number, C, which 
is several orders of magnitude larger than the components of A. A general rule of 
thumb for choosing the value of C is to use C = max\Aij\ x 10 4 . The constant C is 
added to each of the diagonal elements in A^e corresponding to the degrees of freedom 
on the Dirichlet boundary, and C is multiplied by the specified temperature on the 
boundary and added to the corresponding components in the load vector, F^e. For 
example, to enforce the Dirichlet boundary condition on the left end of an element 

A^e = Aqs+C for i = 0, . .. 

Z% 1% 

Fqz = Fqz | Cu s for i = 0 (2.40) 

The Neumann boundary conditions, or natural boundary conditions, in (2.6) are 
included in the right hand side of the weak form (2.10) and are thus included in F n e. 
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Enforcement of the boundary conditions along T D causes the system of equations to 
be nonsingular so that a solution to the global system of equations exists. 

The numerical accuracy of the solution of the global system of equations depends 
on the conditioning of the A matrix and the choice of the value of C. The condition 
number of a matrix is defined as the ratio of the largest eigenvalue to the smallest 
eigenvalue. A matrix is singular if its condition number is infinite since at least one 
eigenvalue is zero, and although the A matrix is not singular due to the enforcement 
of the Dirichlet boundary conditions, a large condition number limits the accuracy 
to which a solution can be obtained. If the value of C in the penalty approach is 
chosen to be large enough, then u « u s on T^. Generally this boundary condition 
is satisfied to machine numerical precision with the proper choice for C and a low 
condition number. Gauss elimination can be used to solve the system of equations 
which involves row reduction to obtain an upper triangular matrix. The system 
matrices are sparse, having a relatively small number of nonzero entries, due to the 
choice of the basis functions. The derivatives of the basis functions are orthogonal to 
each other except for the nodal basis functions, denoted with the 0 and 1 subscript in 
(2.23) and (2.19). The local stiffness matrices, defined as K x and K z in (2.34), use 
the derivatives of the basis functions, and therefore are almost diagonal. There are 
off-diagonal nonzero entries in the first and second rows and columns since the nodal 
basis functions do not have the property of orthogonal derivatives. The zero entries 
in the stiffness and the small number of zero entries in the local mass matrices create 
a sparse A^e matrix when the outer tensor product, defined in (2.31), is used. 

The relatively small number of nonzero entries makes it inefficient to store the 
zero entries. To save memory and time, the global system is solved using a sparse 
gauss matrix solver [21]. 
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Chapter 3 

A Posteriori Error Estimation 


There are essentially two types of error estimates, a priori and a posteriori A priori 
estimates are theoretical bounds on the error in an approximate solution, and establish 
the rate of convergence of the error for a particular method. These bounds are a 
function of the exact solution and the parameters that influence the accuracy of the 
approximate solution. For p- version finite element methods, these parameters are the 
mesh size, /i, and the degree, p, of the polynomial approximation in an element. A 
posteriori error estimates use the approximate solution itself to estimate its error. 
They fall into two categories: implicit and explicit. Explicit error estimates are 
fast and computationally inexpensive as they typically involve post-processing of the 
finite element solution. Implicit error estimates are obtained by solving an additional 
boundary value problem for the error with the approximate solution as input data. 
Implicit error estimates are more expensive than explicit error estimates, but they are 
also more reliable and more accurate. This research focuses on the element residual 
error method, which is an implicit type a posteriori error estimate. The global residual 
is defined as 

R(v) = L{v) — B(u , v) 

It provides information on how close the approximate solution is to satisfying the 
boundary value problem. This chapter discusses the method and approach to pre- 
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dieting the global error in the finite element solution to the two-dimensional boundary 
value problem in (2.12). The global error estimate is a sum of element error estimates 
obtained by solving an element boundary value problem with the element residual 
analogous to a heat source term. The element residual error method developed by 
Ainsworth and Oden [1] for the total error is first described. The performance of 
the error estimate is then demonstrated on two example problems. When using the 
hierarchical modelling approach, the mesh size h, the through-thickness polynomial 
order p z , and the in- plane polynomial order p x are parameters that affect the accuracy 
of the solution. A directional error estimate is proposed to distinguish the modelling 
error (p z ) from the finite element error (p x and h). Its performance is investigated on 
two example problems. 

3.1 Element Residual Method 

The error in the finite element solution is the difference between the exact and the 
finite element solutions, e = u — u. The weak form of the boundary value problem 
for the error is obtained by subtracting B(u,v) from both sides of (2.12). 

B(u,v) — B(u,v) = L(v) — B(u, v) 

B(e,v) = L(v) — B(u,v) 

B{e,v) = R{v) \/veV{Q) (3.1) 

The solution of this global equation yields the actual error, but the solution space 
i/(Q), defined as a subspace of if 1 (Q), is of infinite dimension. The finite element 
approximation described in Chapter 2 for the temperature can also be used to approx- 
imate the solution to (3.1) for the error. However, if the subspace used to approximate 
the error is the same subspace F(fi) that was used for the finite element approxima- 
tion, the result will be that the error is zero. Since the right hand side of (3.1) with 
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v e V^fi) is precisely what was solved for the finite element solution, R(v) = 0 for 
every v e Therefore to obtain a nontrivial approximation for the error, a larger 

subspace must be used. A simple approach to constructing a larger subspace is to 
increase the polynomial degree of the approximation on each element in both the 
through-thickness and in-plane directions. Solving the global error problem on this 
enriched polynomial space would, however, require more computation time than orig- 
inally required to solve the heat conduction problem. A more efficient approach is to 
solve a local (or element) boundary value problem to provide an estimate of the error 
on each element. The local boundary value problem for the error in each element is 
derived from the original differential equation for the temperature. The global differ- 
ential equation (2.6) holds true for any point in the solution domain. Therefore this 
differential equation also holds true for any collection of points defining an arbitrary 
subdomain or element, making (2.6) valid for any element. With the finite element 
solution, ft, in hand, the quantity V T (/^Vft) is added to both sides of (2.6) defined 
on each element 


+V r (^Vft) — V t (k;Vu) 


Q + V t (kVu) 


-V T (nVe) 


Q + V 7 (jvVft) on element Ct e (3.2) 


Appropriate boundary conditions for this equation will be considered with the weak 
formulation, obtained by multiplying (3.2) by a test function v and integrating by 
parts 


I V T (KVe)v dQ e = [ QvdQ e + f V T (nVu)v dtt e 

JQ e JQ e JQ e 


— / (KVe) T n^eV ds — ( 

.JdQ e Jn e 


(nVeVVv dn e 


/ Qv dft e + / (K'Vu) T riQeV ds 
Jn e JdQ e 

[ C KVufVv dn e (3.3) 
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Substituting u = u — e into the boundary integral terms simplifies (3.3) to 

[ (nVefVv dn e = [ Qv dft e - [ (nVufVv dfl e + [ ds 

Jn e Jn e Jn e Jan e 

Bne(e,v) = Fqb(v) - Bne(u,v)+ [ (nV u) T n^v ds (3.4) 

Jdn e 

where F^e (v) is the linear function associated with the source term defined as f Qe Qv dCt e 
The subscript f2 e on B^e and F^e denotes a local domain. Equation (3.4) is the weak 
form of the local error problem. This equation is the element contribution to the 
global error problem in (3.1), since the exact boundary flux terms in (3.4) sum to 
zero along interior element boundaries. The actual error satisfies (3.4), but the exact 
flux in the last term, (/^Vw) r %e, is not known. Therefore an approximation must 
be made for the flux boundary conditions on each element: 

q n e r; -(nVu) T nne (3.5) 

Using this approximation for the exact boundary flux in (3.4) results in the local 
problem to be solved in each element. Find e E V(ft e ) such that 

Bne(e,v) = F Q e{v) - B n e{u,v) - i q^v ds (3.6) 

Jdn e 

for all admissible v E U(f2 e ). The next section discusses two approaches to construct- 
ing an approximate boundary flux, (foe, using the finite element solution. 

3.2 Element Boundary Conditions 

The boundary flux (#vV«) T nne must be determined for each type of element bound- 
ary. Three possible types of element boundaries are considered: Dirichlet boundary, 
Neumann boundary, or interior element boundary. The boundary conditions for the 
edge of an element that is coincident with a global boundary is straightforward. For 
an element boundary that is coincident with T^, where the temperature is specified 
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as a constant, the finite element solution satisfies the essential boundary condition 
exactly 

u = u = u s on dfl e n r D , . 

e = u — u = 0 on d£l e n T D ' ' ' 

and the test function in (3.6) is required to be zero on dfl e n F D . For an element 
boundary that is coincident with T N , the applied heat flux is specified in the problem 
statement as q s . For this natural boundary condition, 

— (KVu) T rtQe = q s on dtt e n T N (3.8) 

The third type of boundary that must be considered is the interior element bound- 
ary (the interface between two neighboring elements) where an approximation to the 
exact flux must be made. The choice of the approximate boundary flux determines 
the accuracy of the error estimate. The next sections discuss two methods of approx- 
imating the exact flux on the interior element boundaries, the average flux method 
and the equilibrated flux method. 

3.2.1 Average Element Boundary Flux 

The finite element approximation results in a flux that is discontinuous across inte- 
rior element boundaries. One approach to approximating the exact flux on an interior 
element boundary is to average the fluxes computed from the finite element approxi- 
mation on the boundary between neighboring elements. Let Q e/ denote the neighbor 
to element and let 7 denote the shared edge, then the average flux on the shared 
edge is given by 

1 j 7 

?n«| 7 = 2 (-«V«ne| 7 - «V«ne/| 7 ) n Q e (3.9) 

where n^e is the outward unit normal vector to the edge of element fi e . The aver- 
age flux, applied to the shared boundary of the neighboring elements, satisfies the 
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continuity condition 


<7n e l 7 + <7o e, | 7 = 0 on the shared edge 7 = dd e D dQ el 
since = —n^ and 

- 1 ( x A \T 

?fi«'l 7 = - (-«V%'| 7 - «Vtt n *| 7 ) n Q e> 

1 / , \T 

= -- (-KVV | 7 - «VttQe| 7 j njje 

= - tfn <=| 7 

3.2.2 Equilibrated Flux 


(3.10) 


(3-11) 


Another method for approximating the fluxes on the interior element boundary is to 
develop boundary fluxes that, in addition to satisfying continuity (3.10), satisfy an 
equilibrium condition for that element [1]. The equilibrated flux method as described 
in [1] is applicable to general two-dimensional meshes of quadrilateral and triangular 
p-elements. It simplifies considerably for hierarchical modelling. 

The equilibrium condition for an element is derived by integrating (2.6) over an 
element 


[ - V T (#eVu) dn e = [ Q dQ e 
Jn e JQ, e 

[ —(K'Vu) T riQe ds = [ Q dVt e 

JdQ e JQ e 


(3.12) 


In the second step, integration by parts is used on the left hand side. For steady state 
conditions, the flux around the boundary of an element must be in equilibrium with 
the internal heat source. Therefore, the approximate flux must be constructed such 
that for every element in the mesh 



Note that (3.13) is equivalent to setting v = 1 in (3.6), that is, 

0 = Fne(l) - B Q e(u, 1) - [ q Q e ds 

JdQ e 


(3.13) 


(3.14) 
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since Ffoe(e, 1) and Ffo e (u, 1) are zero. The term Bne(u , 1) is retained in (3.14) because 
it is part of the element residual and plays an important role in constructing the 
equilibrated approximate flux. The equilibrium condition (3.14) can be extended 
to higher orders by replacing 1 in (3.14) with a set of two-dimensional nodal basis 
functions, 0 n , defined at each node n. 

Fn*(O n ) - B n e(u,d n ) - / q^e n ds = 0 forn=l,...,n 7 (3.15) 

JdQ e 

The On s represent basis functions associated with the nodes on the boundary of an 
element. The specific form and the required number, n 7 , of the nodal basis functions 
along an edge, to be discussed later, will depend on the approximation selected for 

fa*- 

The approximate flux must satisfy the equilibrium condition (3.15) along with 
the continuity requirement (3.10) and the Neumann boundary condition (3.8). For 
hierarchical modelling of a two-dimensional problem, there are only one-dimensional 
elements in the ^-direction. In order to apply the equilibrated flux method to the 
hierarchical model, an element is conceptually expanded in the through-thickness 
direction into a rectangular element. The node and edge numbers are defined on 
the expanded local element as shown in Figure 3.1. For this expanded element, the 
top and bottom boundaries of each element are coincident with the global domain 
boundary, T, and therefore an approximate flux is only required for the left and right 
boundaries, denoted by 72 and 74 in Figure 3.1. 

The approximation for the left and right element boundary flux is assumed to 
be a polynomial of degree p z , the same as the polynomial degree of the hierarchical 
model. The approximate flux on an edge, 7, of the boundary of element is written 
in terms of the nodal basis functions, 0, as 

Pz + 1 

(3.16) 

i— 1 
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Figure 3 . 1 : Convention for node and edge labels for an expanded local element. 


where i represents the degrees of freedom on the edge 7 and a* are unknown coeffi- 
cients. Due to this choice for the approximate boundary flux, the number of nodal 
basis functions required for the left edge or the right edge is n 7 =p z +l. The complete 
set of two-dimensional basis functions, 0 , is given by 

MOMv) } n ode 1 

MOMv) } node 2 

MOMv) } node 3 

MOMv) } node 4 

MOMv) 

. > edge 71 

MOMv) 

MOMv) , 

0 = . f edge 72 ( 3 . 17 ) 

MOMv) 
foiOMv) , 

. > edge 7s 

<f> Px (OMv) > 

MOMv) 

MOMv) 

. > edge 74 

. MQ<t> P Sv) > 

The two-dimensional basis functions, are based on the finite element basis functions 
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used for the finite element approximation in the ^-direction (2.23). They are essen- 
tially a tensor product of one-dimensional basis functions for a quadrilateral element. 

There are p z + 1 two-dimensional basis functions associated with the nodes of an 
edge, which is where the flux is being approximated: one for each vertex node and 
p z — 1 for the edge node. Therefore there are p z + 1 equations that must be solved 
to determine the coefficients ct; in (3.16). The set of two-dimensional basis functions 
associated with the nodes of edge 7 is a subset of the complete set given in (3.17 and 
defined as 




'jp I vertex nodes of 7 

e\ ' 

ej 

> edge nodes of 7 


Q 1 

U Pz + 1 


(3.18) 


The first four two-dimensional basis functions for edge 72 are shown in Figure 3.2. 



Figure 3.2: First four two-dimensional basis functions for edge 72. 


The last term in (3.15) is an integral over the entire boundary of an element, and 
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can be written as a sum of integrals over each edge of the element: 


Fae(6 n ) — Bae(u,O n ) - 


[ Qb^n 1 ds + f Qqc ds + f q t 6 7 ds 

_J 7 1 J 1 2 ^73 


+ [ Qn e 6Z 4 ds 

^74 


73 

= 0 


(3.19) 


where q t is the heat flux specified on the top boundary of an element, and is the heat 
flux specified on the bottom boundary of an element since 71 and 73 are coincident 
with T n . Substitution of (3.16) into the boundary integral terms in (3.19) which 
include the approximate boundary flux, yields 


J^ei ds = J yYl a i e ij 


&Z ds for n = 1 , . . . , p z + 1 


(3.20) 


This equation is true for each node, n, along 7 , and can be rewritten in terms of a 
mass matrix as 


/ q^eO 7 ds = M 7 a (3.21) 

J 1 

where M 1 is the mass matrix along edge 7 with entries 

Mij = [ QiOj ds for ij = 1, . . . ,p z + 1 (3.22) 

In order to solve for the coefficients, a, (3.19) is written using (3.21) and moving the 
other terms to the right hand side. Due to the properties of the nodal basis functions, 
the sum of boundary integrals in (3.19) is simplified. For a nodal basis function at a 
vertex node, there are two edges along which the basis function is nonzero. Whereas 
for a nodal basis function at an edge node, there is only one edge along which the 
basis function is nonzero. This can be seen in Figure 3.2 where the first two plots 
show the basis functions for 72 , and the third and fourth plots show the edge basis 
functions for 72 . Each of the higher-order basis functions associated with the edge 
node is nonzero only on the edge containing the edge node. Using this information 
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about the nodal basis functions, and substituting (3.21) into (3.19) results in the 
following equations. For 72, 


For 74, 


M 12 a 


' F Q e(d f) - B n e(u, 6?) - / 7l q b 0? ds ' 

F^{6f) - B Q e(u,d?) - f 73 q t e? ds 
F,Aef)~B^(u } 0f) 1 


(0, 


172 

Pz + lJ 


-B< 




Q12 > 

U Pz + 1J 


M 74 ct 


' F n e(d?) - B Q e(u, or) ~ f 71 q h 0 T ds > 
Fn*(e?) - B n *{u, Of) - f J3 q t 0f ds 
Fne(0?)-B Q e(u,9?) 


(3.23) 


(3.24) 


{ F^e^-Bne^e;^) J 

Equations (3.23) and (3.24) are used to determine the approximate flux on each 
respective edge. The system of equations that must be solved to determine the coef- 


ficients a will always have a unique solution since the mass matrix in these equations 
is positive definite. 


3.3 Error Indicator 


The error residual problem, shown in (3.4), can be used to approximate the error 
after choosing the method of approximating the true flux on the element boundaries 
and choosing an appropriate subspace V(Q e ). The approach to approximating the 
error is the same as the approach used for the finite element approximation in Section 
2.3.3. The approximation of the error in an element is 

Px+G Pz+<T 

e = Y, E MxWjWbn (3.25) 

i=0 j = 0 

The Bubnov- Galer kin method is used so that the test functions, v, are the same as 
the approximating functions, 

v = Xi i = 1, • • • , ip x + 1 + cr) x {jp z + 1 + <r) (3.26) 
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Upon the substitution of (3.25) and (3.26) into (3.4), the approximate solution to the 
error residual problem for each element becomes a system of equations and can be 
written in matrix form as 

Bne(e,v) = Fne(v) — Bne(u,v)+ [ qnei) ds (3.27) 

JdQ e 

A e b = F e (3.28) 

where A e represents the coefficient matrix based on the chosen subspace, V^c, b is 
a vector of unknown constants, and F e is the load vector consisting of the element 
residual term and the boundary flux term. This system of equations is solved for the 
unknown coefficients, b. The finite element solution satisfies the essential boundary 
condition on T^. Therefore the error must be zero on this boundary. The penalty 
method, described in Chapter 2, is used to enforce this essential boundary condition 
in (3.28) for elements with boundaries coincident with Tp. This system of equations 
is singular except on those elements with a Dirichlet boundary where the error is set 
to zero. Singular value decomposition is used to find the the solution for b which 
minimizes r = \A b—F\ where the number r is called the residual of the solution [21]. 
In other words, the least-squares best compromise solution is found. The singular 
values are zeroed out when they are less than some tolerance. Therefore discretion 
in choosing a value for the tolerance must be used. 

The error indicator for each element is obtained by calculating the energy norm 
of the error function. The energy norm is defined by the weak formulation 

£n e = IPHIne (3.29) 

i 

= [ («Ve) T Verffi 2 (3.30) 

J Q 

The global error estimate is calculated by computing the sum of the squares of the 


'Bne (e,e) 
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local error indicators 


£ = 


n 1/2 


£(^) 2 

_Q e 


(3.31) 


3.4 Numerical Examples for the Total Error 


The performance of the error estimates will be examined in this section. Example 
problems are chosen such that the exact solution is known. The problems are not 
meant to be realistic problems. The use of problems with exact solutions allows 
the actual error in the finite element solution to be computed and compared to the 
estimated error computed from (3.31). The first example problem has a smooth 
solution composed of a sine and a cosine function. The second example problem has a 
rough solution composed of an arctangent function. The finite element approximation 
for the first example will generally provide a better approximation (less error) than the 
second example problem for the same orders of approximation and mesh size because 
of the smoothness of the exact solution to the first example problem. The problems 
are developed using an inverse approach by first choosing a closed form solution, u, 
which satisfies the desired boundary conditions. Then the internal heat generation, 
Q, required to obtain the closed form solution is determined by substituting the exact 
solution into the heat conduction equation (2.6). 

The examples will consider the two dimensional heat conduction problem shown 
in Figure 3.3. The heated plate has dimensions L by d with the following boundary 
conditions: 


u 

du 

' dz 
du 

' dz 


0 along the edges x = 0 and x = L 

0 along z=— ^ 

£ 

qt along z = ^ 


(3.32) 


The plate has an internal heat source 0. and the heat load q t applied along the 
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boundary at z = f is a function of x. The bottom of the plate is insulated, and 
the sides are held at a constant temperature of zero. For all results presented in this 
section, k x = k z = 1.4, L = 1, d = 1, and the value of a = 1 in (3.25) are used. 



Figure 3.3: Two-dimensional plate for numerical examples. 


3.4.1 Case I: Smooth Solution 


The exact solution to the heat conduction equation (2.6) is chosen to be the analytic 
function 


u 


a 



—2d . (i tz 

- — sin — 

Kz'K \ 2d 



(3.33) 


This closed form solution satisfies the boundary conditions in (3.32) with the following 
heat load applied to the top boundary 


q t (x ) = 1 — cos 



(3.34) 


The required internal heat source determined from (2.10) is 


Q(x,z) 


cnr k z L 2 — (16d 2 k x + k z L 2 ) cos(^) sin[ ^ rf 4rf 2 ^ ] 
2 dkj? 


(3.35) 
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The contours of the exact solution for the heated square plate with a = 10 is shown 
in Figure 3.4. The solution is symmetric about the center line at x = 0.5 and exhibits 
a very smooth behavior. 



Figure 3.4: Exact solution for Case I. 

The error estimate is computed using the average flux approximation, (3.9), and 
compared with the actual error in the energy norm in Figure 3.5. The graphs in Figure 
3.5 show the rate of convergence of the error in the finite element approximation with 
respect to the length of the elements, /i, in a uniform mesh. The error in the energy 
norm is plotted versus the inverse of the mesh size, which for a uniform mesh is the 
same as the number of elements, so that the error curves have a negative slope. Each 
graph represents a fixed particular value for p x . Each line in a graph represents a 
different value of p z varying from one to six. Each point on a particular line represents 
a uniform mesh of ^ one-dimensional elements. As p x increases, the error converges 
faster (with fewer elements). The rate of convergence of the error, indicated by the 
slope of the solid line in each plot is equal to the value of p x . For the constant part 
of a convergence curve where the lines level off, the hierarchical modelling error is 
dominating. In this region, refinement of the mesh does not improve the finite element 
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solution, and an increase in the chosen value for p z is required to improve accuracy. 
When the hierarchical modelling error dominates, the estimated error (denoted by 
dashed lines in Figure 3.5) is only slightly larger than the actual error for all values of 
p x . When the finite element error dominates, difference between the estimated error 
is significantly larger than the actual error, but exhibits correct rates of convergence. 
The discrepancy between the estimated and actual error is more pronounced for even- 
orders of p x . 

The error estimate computed using the equilibrated flux is compared to the actual 
error in Figure 3.6. Comparing Figure 3.6 with Figure 3.5 shows that the equilibrated 
flux method provides more accurate error estimates, particularly when the error is 
dominated by the finite element approximation error. The estimated and actual errors 
using the equilibrated flux method are virtually indistinguishable for all values ofp*, 
p z , and mesh size. 


38 






estimate 


Figure 3.6: Error for Case I using the equilibrated flux method. 
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3.4.2 Case II: Rough Solution 


The performance of the error estimates for a problem with a rough solution will now 
be examined. The closed form solution is given by 


u = a(L — x ) [arctan [a{x — rr 0 )] + arctan (aaro)] 


2d . (i xz 7 r 

sin 


k z 7T 


2d 


(3.36) 


This solution satisfies the boundary conditions (3.32) with the applied heat load along 
the top of the plate defined by 


qt(x) = a(L — x ) [arctan [ a(x — Xq)] + arctan(c>!Xo)] (3.37) 


The internal heat source required to deliver the solution in (3.36) is determined from 
(2.6) to be 


Q(x,z) 


aiv(L — x) [arctan[a(a; — Xo)] + arctan[<xEo]] sin 

/ 7 T 7 TZ \ 

\ 4 2d) 

2 adk x 

2d 

2asin(f — ff ) 2a 3 (L-x)(x-x 0 ) sin(f--ff ) 


7t[1+q; 2 (®— *q) 2 ] 7t[1+q; 2 (*-*o) 2 ] 2 


k z 


(3.38) 


where a = 30 and £ 0 = 0.3. This exact solution has steep gradients in the ^-direction 
in the vicinity of x = Xo as seen in the contours of the exact solution shown in Figure 
3.7. 

The actual error is compared with the estimated errors using the average and 
equilibrated boundary flux methods in Figures 3.8 and 3.9. The actual error in the 
energy norm is not as well behaved as in Case I because the solution for this case 
is not as smooth. Smaller mesh sizes are required to reach the expected rate of 
convergence of p x [4] as shown in Figures 3.5 and 3.6. The error estimates obtained 
using the average flux method are much larger than the actual error. In comparison 
of the error estimates obtained using the average and equilibrated flux methods for 
Case II, the equilibrated flux method provides noticeably better error estimates than 
the average flux method. A contour plot of the pointwise error in a finite element 
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Figure 3.7: Exact solution for Case II. 


solution obtained with 32 elements and p x = p z = 2, shown in Figure 3.10, shows 
that the greatest errors occur after x = Xq where the steep gradients occur. The 


element error indicators, obtained with the average and equilibrated flux methods 
are compared with the actual element error in Figures 3.11 and 3.12. The actual 
and estimated errors in the finite element solution for p x = p z = 2 with 32 elements 


are shown. For this case, the global error for the equilibrated flux provides a much 
better estimate of the error than the average flux method. The average flux method 
severely over-estimates the error for each element in the region near x = Xo, while the 


equilibrated flux method provides an accurate indicator of the error for each element. 
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estimate 

Figure 3.8: Error for Case II using the average flux method. 
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actual 


estimate 

Figure 3.9: Error for Case II using the equilibrated flux method. 
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Figure 3.10: Pointwise error in the finite element solution for Case II using 32 elements 
with p x =p z = 2. 


890 0 




0 0.25 0.5 0.75 1 

X 


Figure 3.11: Element error indicators for Case II solution using 32 elements with 
p x = Pz = 2 using the average flux method. 
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Figure 3.12: Element error indicators for Case II solution using 32 elements with 
p x = Pz = 2 using the equilibrated flux method. 
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3.5 Directional Error Indicators 


The element error indicator described in the previous section gives an estimate of the 
accuracy of the approximate solution on each element. There are various approaches 
to improve the finite element solution. Some of these techniques include refining the 
mesh, increasing the order of the approximation, and using different basis functions. 
When using the hierarchical modelling approach, the through- thickness polynomial 
degree, the mesh size, and the in-plane (or x-direction) polynomial degree are pa- 
rameters which will affect the accuracy of the solution. A directional error indicator 
provides an estimate of the error in the hierarchical model and finite element ap- 
proximation separately. The directional error indicator indicates which parameters 
to change to efficiently improve the approximate solution. According to Figure 3.6, 
the error in approximating the solution in the ^-direction with a polynomial of order 
p z is dominant at the points where the lines are fiat. In this circumstance, increasing 
the value of p z will improve the approximate solution while refinement of the mesh 
will not. If the approximate solution has an error which is located on the part of the 
lines where the slope is p x , then the error in the finite element method is dominant, 
and refinement of the mesh would be the appropriate choice to reduce the error in 
the solution. Directional error indicators could be used to distinguish between these 
circumstances by determining which contribution to the error, the hierarchical model 
or the finite element method, is the largest. 

The error in the solution is composed of the error in the hierarchical model, e Model, 
and the error in the finite element approximation of the hierarchical model, 6fe- 

e Model = U — U Model (3.39) 

6fE = U Model — U (3.40) 

The solution UModei is the exact solution to the hierarchical model, that is the exact 


47 



solution to the system of differential equations resulting from the assumption that the 
solution is a polynomial of p z in the through-thickness direction of the entire domain. 
The finite element error is associated with the polynomial approximation pf degree 
p x in the x-direction. The total error can be expressed as the sum of the modelling 
error and the finite element approximation error 

e = u — u 

& 'Ll 'LL Model T 'LL Model LI 

£ = 6 Model + &FE (3-41) 

The energy norm, defined in (3.30) from the weak formulation, is used to measure 

the error in the approximate solution 

|||e||| 2 = B(e, e) (3.42) 

Substituting the definition of e in (3.41) into (3.42) and expanding yields 

B(e, e) = B(e Model + &fe, e Model + £fe) 

B(e,e) = B(e Model, e Model) + B(epE , £Fe) + B (e Model, &fe) + 

B(eFE , e Model ) (3.43) 

i?(*, ♦) is symmetric, and the modelling error is orthogonal to the finite element ap- 
proximation error [5] , meaning 


B(c Model, &Fe) ~ 0 

(3.44) 

B(&FE, ^Model) = 0 

(3.45) 


Using (3.44) and (3.45), the last two terms in the definition of the energy norm of 
the total error in (3.43) become zero, so that 

B{e,e) = B{e Model, e Model) + B{6 fE, £fe) 

1 1 Mil 2 = 1 1 \ & Model 1 1 P T || 1 1 2 (3.46) 
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This equation shows how the directional errors are related to the total error. 

The total element error for the model problem is approximated using the following 
polynomial, as discussed in Section 3.3, 

e = E E Mx)^j(z)bij cr > 1 (3.47) 

i = 0 ,7=0 

Recall that a must be greater or equal to 1 to ensure that the result for the error 
estimate is nontrivial. The finite element error is the error in the solution which is 
controlled by the mesh size and the polynomial order in the x-direction, h and p x . 
The modelling error is controlled by the through-thickness polynomial order, p z . To 
approximate each of the directional error indicators, the opposing a is set to zero in 
(3.47). The directional error polynomial approximations are 

Vx+V Pz 

&FE = X X 
i=0 7=0 

Px Pz+<J 

& Model = EE <j)i{x)^j(z)Cij (3.48) 

i=0 j=0 

Using each of these approximations separately in the element residual error equation 
(3.28) with either the average and equilibrated flux methods described in Sections 
3.2.1 and 3.2.2, the approximate solution for each of the errors can be found. The 
directional error indicators are computed using the energy norm of the approximate 
error functions: 


2 

n 


£FE 


£ Model 


B(e,e) = ^#ne(e,e) 

Q e 

|||ei?£;|||n = B(sfe, Gfe) = . ;E Bge ( cfb, &fe) 

V 

1 1 1 GjVf orfei 1 1 |fi = \j B(e Mode i,e Model) = , ;E Bge (e Model , & Model ) 

V Q e 


(3.49) 

(3.50) 

(3.51) 
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3.6 Numerical Examples for the Directional Error 
Estimates 

Results for Case I and II are examined using the directional error indicators. The 
equilibrated flux method is used to approximate the element boundary fluxes. Figure 
3.13 shows the actual error, the model error estimate, and the finite element error 
estimate plotted versus the inverse of the mesh size for the problem in Case I. The 
p z value is varied with uniform mesh refinements for a fixed p x value. The graphs 
in Figure 3.13 show that the estimated error in the finite element approximation 
(denoted by triangles) has a slope of and follows the convergence rate of the actual 
error. The actual hierarchical modelling error for each value of p z is determined by the 
horizontal portion of the actual error line. The estimated hierarchical modelling error 
is accurate for sufficiently small mesh sizes, which are determined by the value of p x . 
At each point on each line, the most effective method of improving the approximate 
solution is determined by which error indicator is larger. 

The directional error indicators for Case II are shown in Figure 3.14. The actual 
error, model error estimate, and finite element error estimate are plotted versus the 
inverse of the mesh size. The actual error norms are not as well behaved as in 
Case I due to the roughness of the problem, however, the finite element error is 
accurately predicted by the directional error estimates. The estimated hierarchical 
modelling error is slightly higher than the actual hierarchical modelling error for low 
values of p z and for large mesh sizes. The results of these examples show that the 
directional error indicators accurately determined which error dominates, the model 
or the finite element error, for a smooth and a rough solution. Knowledge of which 
error is dominant in the approximate solution indicates the most effective way to 
improve the approximate solution. 
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Figure 3.13: Directional error estimation for Case I using the equilibrated flux 
method. 
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Figure 3.14: Directional error estimation for Case II using the equilibrated flux 
method. 
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Chapter 4 

Extension to Multi-layered 
Materials 


Structural concepts for aerospace vehicles often consist of panels made of multi-layer 
materials. Examples include skins of composite laminates or metallic thermal pro- 
tection systems. Accurate and efficient analysis methods are required to aid in the 
design of such concepts. Hierarchical modelling, p-version finite elements, and error 
estimation are extended to multi-layered materials in this chapter. 

4.1 Boundary Value Problem 

Section 2.1 presents the boundary value problem for a single layer. The two-dimensional 
domain for a multi-layered material for a plate with length L, thickness d, and n /, lay- 
ers of material is shown in Figure 4.1. The boundary value problem must be satisfied 
by the solution for the temperature in each layer, Ui, with the boundary conditions 
described on the boundary of the global domain: 

-V r (K,: Vui ) = Q i = l,...,riL 
-(KiVui) T n n = q s on T N 

Ui = u s on T/j (4.1) 

For layered materials, continuity at the boundaries of the layers must also be satisfied. 
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Figure 4.1: Two-dimensional domain for a multi-layered material. 


At each layer interface, the temperature and the heat flux are required to be equal 

i = l,...,n L ~l (4.2) 


U i 


— (Ki'Vui) 1 n ; 


u i + iU 

— («*+ 1 V Ui+i ) r n i+ 1 1 


where i denotes the layer of interest, z ? - represents the z-coordinate of the layer inter- 
face, and the outward unit normal vector to the boundary of layer i is n^. 


4.2 Finite Element Method for Multi-Layers 


The finite element method used in this research for materials with multiple layers 
is similar to the method described in Chapter 2. The through-thickness direction is 
modelled by a polynomial function of degree p z , reducing the spatial dimension of 
the domain by one, a mesh is created for the dimensionally- reduced domain, and the 
p-version finite element method is used to complete the approximation. The thermal 
conductivity for the orthotropic multi-layered materials is assumed in this research 
to be a piecewise constant function in the through-thickness direction, 


(kx)i 0 
0 (k z )i 


(4.3) 
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4.2.1 Weak Formulation 


The weak form of (4.1) is constructed by multiplying by a test function and integrating 
by parts. The test function is required to be zero on T c . 

f (nVufVv dQ= [ QvdQ- [ q s v ds (4.4) 

Jo, Jq Jr N 

In abstract form, the problem is to find u e such that u = u s on F D and 

B{u , v) = L(v) for all admissible v e V^fi) (4.5) 

The subscript i has been dropped in (4.4) and (4.5) with the understanding that the 
thermal conductivity is a piecewise constant function of z, and the exact solution is 
piecewise defined on each layer. The hierarchical modelling combined with p-version 
finite element method is used to approximate the solution to (4.5). The p-version 
method for multi-layers is identical to the method described in Section 2.3.2 for a 
single layer, where the basis functions, 0, in (2.23) are used. However, since one 
approximate solution is computed for all layers in each element, the hierarchical 
modelling method is modified by choosing a more efficient set of through-thickness 
basis functions to approximate the exact solution using piecewise polynomials. 

4.2.2 Optimal Hierarchical Model 

Piecewise constant thermal conductivities in the boundary value problem in (4.1) 
result in a piecewise exact solution in the through-thickness direction. The work 
of Vogelius and Babuska [5] shows that the optimal choice of basis functions for 
hierarchical modelling of the boundary value problem (4.1) are piecewise polynomials. 
The choice of basis functions are optimal in the sense that the approximate solution 
converges to the exact solution at the expected rate as the polynomial order, p z , 
approaches infinity and as the plate thickness, d, approaches zero. The basis functions 
for hierarchical modelling in single-layered materials are referred to as homogenized 
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basis functions. The set of optimal basis functions, rj>, is constructed by scaling the 
odd-degree through- thickness homogenized basis functions, Vb by the conductivity 
of each layer and shifting the resulting polynomial to enforce continuity at the layer 
interfaces: 


where 


to 

II 

to 

05. 

0 — 
2 


4-1 

— 1 < V < Vi 

to 

1 

4-i 

Vi <V <f}2 


. 4^ 

Vn L - 1 <!/<!, 



(4.6) 


= Wh Pv - 1 

4-i = (Vi-i)) + Ci i = 2,...,n L 

C2 = (Ej7^- l( " ,) 

c i+ 1 = Ci + - 1>2j-i(Vi-i)) i = 2, • • • , n L - 1 (4.7) 

\™z )i 

and i is the layer of interest, is the optimal polynomial of degree 2 j — 1 in the 
i th layer, and ipj are given in (2.19). These basis functions provide more accurate 
approximations for layered materials than the single layer through-thickness basis 
functions described in Section 2.3.1. An example of these basis functions are shown 
in Figure 4.2 for a three layer material with (k z ) i = 1, (k z ) 2 = 10, and (k z ) 3 = 1 for 
the polynomials of degree 0 through 7. Note that the odd degree basis functions have 
discontinuous derivatives at the layer interfaces. 


4.2.3 Element Equations 


The finite element solution, u, on an element is written as a linear combination of the 
basis functions 

Px Pz 

= = fp) T a e = X r a e ( 4 . 8 ) 

i=0 j = 0 
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Afte = 


Vx(kVx) t dQ e 


JQ e 

n L 


dzdx 


,v l> nx\ nzi 

E / / Vx(kVx) 5 

i=1 7*0 J Zi_\ 

fZi ^ 

dx ®22 ( k x )iipip dz + 

i = 1 


r x i d<p d<f> 

lx o <9x <9a; 

rx i 




, ,T , ^ST' f Zi n \ dtp dtp 

*+ ix °Z L^WaT 


dz 


Fne = 


[ Qx d£l e - [ q s X ds 
Jn e Jdn e nr N 

f*x i rzi /* 

E / / dzdx - / ds 

i=i *' x o J Zi~ i 7a^ e nrjv 


(4.11) 


(4.12) 


4.3 Multi-Layer A Posteriori Error Estimation 


The element residual method is applied to multi-layered materials. A local boundary 
value problem for the error is developed, and the equilibrated flux method is mod- 
ified for application to multi-layered materials. The solution to the error equation 
is approximated and used in the definition of the energy norm to compute the error 
indicator for each element. Numerical results are presented for multi-layered materi- 
als. According to the literature search performed, numerical results for applying the 
element residual method to multi-layered materials have not yet been presented. 

4.3.1 Element Residual Method 


The weak form of the local error residual problem derived in Section 3.1 is used for 
the multi-layer problem. 


[ (AVe) 
Jn e 


T VvdQ e = 


[ Qv dQ e — [ (nVufVv dtt e + [ (K'Vu) T nnev ds 
JQ e 7Q e JdQ e 


B Q e(e,v) = Fne(v) - Bne(u,v) + (KVii^rineV ds 


e = 0 onr c (4.13) 

The last term in the variational form of the element error residual problem is the exact 
flux, and an approximation must be made since the exact solution is not known. The 
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next section discusses the modified equilibrated flux method of approximating the 
boundary fluxes on an interior element boundary for multi-layer problems. 

4.3.2 Equilibrated Boundary Flux Approximation 

The equilibrated flux method is used to approximate the element boundary flux, but 
the method must be modified for multi-layered materials. The exact flux in the re- 
direction is discontinuous along the boundary of an element (not across the boundary) 
in the through-thickness direction because of the piecewise constant thermal conduc- 
tivity. Figure 4.3 shows an example of the discontinuous boundary flux on an element 
with four layers. Therefore the method chosen to approximate the exact boundary 
flux should also be discontinuous. To develop an equilibrated discontinuous flux, the 



Figure 4.3: Example of a discontinuous exact boundary flux in the x-direction on a 
four layer element. 

equilibrated flux method is modified by defining the approximate flux on each element 
as 

Qn e = Qn e + Qn e (4-14) 
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where is the average flux and q^e is the deviation from the average required 
to equilibrate the flux. The average flux is discontinuous since it is computed by 
averaging the finite element fluxes from neighboring elements as described in Section 
3.2.1 

1 T 

qn ‘ = 2 OV&Qel^e + KVfinvIan.) n n e (4.15) 

The equilibrium condition for the equilibrated flux in a single layer is used to find 
qne. The definition for q^e in (4.14) is substituted into (3.15) 

Fne(O n ) - B Q e(u,e n ) - / qoeO n ds- qo.e6 n ds = 0 

Jd O e Jd Q« 

for n = 1, . . . ,n 7 (4.16) 

For the boundary of an element that coincides with Tjv, 


qne + q n , = q s on dQ e n Tjv 


(4.17) 


The approach to solve (4.16) for is the same as described in Section 3.2.2 to solve 
(3.15). The deviation from the average flux, ^e, is written in terms of the nodal basis 
functions, 0 n , (see (3.18)) with the same through-thickness polynomial degree used 
in the hierarchical modelling of the plate, p z , as 


Pz + 1 

^?o e l 7 ^ ^ ^ 4^(7 (4*18) 

i= 1 

Equation (4.18) is substituted for q^e in (4.16), resulting in a system of equations for 
the unknown coefficients. For the edge 72 , 


M 12 cx 


Fne(0?) - B Q e(u, Of) - f 72 q n e0T ds - f 71 q b 0f ds ' 
F Q e(0?) - B n e(u, Of ) - / 72 q Q e6 ] 2 ds - / 73 q t 0f ds 
Fne(Or) - B Q e(v, Of) - f j2 qneOf ds 

, Fn-K+i) - B n *(u,6;i + i) - / 72 &< + i ds J 


(4.19) 
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For 74 


' F Q e(d^) - B Q e(u,6 ?) - / 74 q^OT ds - / 71 q b 0T ds ' 

F,A6f ) - B^{u,6f) ~ / 74 ds - / 73 q t 6f ds 

M 74 « = \ F n .{6f) - F>Qe(u, 6f) - / 74 qn*6? ds (4.20) 

, FcAOlUi) - Bw(u,ei + i) - 4 ^<+1 ds , 
where M 7i is a mass matrix along 7 ; defined in (3.22). Equations (4.19) and (4.20) 
are solved for n: ; . which determines the deviation, <jp e - The solution for f/o« and the 
definition for the average flux in (4.15) are substituted into (4.14) to determine the 
equilibrated flux applied to the interior element boundaries. 

4.3.3 Multi-Layer Error Indicator 

The error function in the error residual problem (4.13) is approximated by using the 
equilibrated flux method for in the last term. The approximate error function, e, 
is defined as a linear combination of the basis functions used in the finite element 
approximation 

Px+crpz+cr 

e=Yl (4.21) 

i=0 j = 0 

where a > 1. The definition for e is substituted for e and (3.26) is substituted for v 
in (4.13), creating a system of equations that is written in matrix form: 

Bw (e, ft) = F Qe ( v ) - B Qe (ft, v) - q^v ds 

JdQ e 

A e b = F e (4.22) 

The error indicator for each element is determined using the energy norm of the 
approximate error function e defined as 

= IPIHoe = Vi?ne(e,e) = [ [ (/tVefVe do\ 2 (4.23) 
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The global error indicator is calculated by computing the sum of the squares of the 
local error indicators 


£ = 


-| 1/2 


E (^) 2 


4.4 Numerical Examples 


(4.24) 


Results are shown for a two- layer plate with length L = 2, thickness d 1 , no internal 
heat generation, ( k x )i = (k z ) i = 1, (k x )2 = (k z ) 2 = 10, and a = 1, shown in Figure 
4.4. The layers are of equal thickness, and the boundary conditions are 


u 0 along the edges x = 0 and x = L 


, du 
-kz-z~ 

oz 


-k: 


du 

dz 


0 along 2 = 


4q 0 x(L — x) 

V- 


_d 

”2 

along z 


d 

2 


(4.25) 


where qo = —10. The exact solution is an infinite Fourier series defined for each layer 



L=2 

Figure 4.4: Two-dimensional plate with two layers for numerical results. 


oo 

W =Esm 

n=l 


nnx\ 

~r ) 


A n sinh {a l n z) + B n cosh (a l n z) 
C n sinh {a 2 n z) + D n cosh (ofcz) 


~i<z<o 

0 < -2 < f 


(4.26) 
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where 


C n 

D n 


n 7r 


Oil = 


A„ = 


B„ = 


{kx)i 


32 Lq 0 tanh (Aa l n 


n 4 7r 4 y^M^)2 [sinh (f a*) + cosh (f o&) ^ tanh (f o&) 
A„ 


tanh (f a* 


— An 


(fcx)l (&; 


•z) 1 


\ {KUK) 2 

^4r>. 


tanh ( |q4 


(4.27) 


The contours of the exact solution for the plate problem is shown in Figure 4.5. 



X 


Figure 4.5: Exact solution for the two layer example. 


To illustrate the improvement in the accuracy of the finite element solution ob- 
tained with the optimal hierarchical model, the actual error obtained using the new 
through-thickness basis functions is compared with the actual error obtained using 
the homogenized basis functions in Figure 4.6. Each graph shows the actual error in 
the energy norm versus the inverse of the mesh size for a fixed value of p x . Each line 
in each graph represents a fixed value of p z with uniform mesh refinement. For each 
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value of p x , the results show that the error is less when using the optimal through- 
thickness basis functions except when p z = 1. The rate of convergence of the error is 
p x for the optimal basis functions as seen in the first graph in Figure 4.6. The hier- 
archical modelling error becomes dominant after only a few mesh refinements when 
using the single-layer (homogenized) basis functions. The hierarchical modelling er- 
ror is much larger for the homogenized basis functions than for the optimal ones for 
the same model order, p z . Also, the hierarchical modelling with the homogenized 
basis functions fails to converge with p z as the modelling error does not decrease 
significantly as p z is increased. 

The performance of the error estimate for the multi-layer example is shown in 
Figure 4.7. The comparison of the actual error and the estimated error shows that 
the error estimate predicts the error in the solution reasonably well. 
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o Optimal hierarchical model 


Figure 4.6: Actual error in solutions obtained for two layer example. 
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actual 

estimate 


Figure 4.7: Comparison of the actual and estimated error for the two layer example. 
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Chapter 5 

Concluding Remarks 


An overview of the research presented in this thesis is discussed in this chapter. A 
summary of the main topics of this paper is first presented followed by the conclusions 
resulting from this research. Finally, recommendations for future work in this area of 
research are given. 

5.1 Summary 

Finite element methods are widely used in engineering analysis of many types of 
problems. Hierarchical modelling combined with p-version finite elements is an ad- 
vanced method of approximating solutions to a number of engineering problems, and 
are applied to heat conduction in this work. Error estimation provides a qualitative 
measure of the accuracy of the finite element solution. A posteriori error estimation 
is necessary for adaptive schemes, where mesh refinement or polynomial enrichment 
is automatically controlled to efficiently improve the accuracy of the solution. For 
hierarchical modelling combined with the p-version finite element method, additional 
information about the error can be determined through the use of directional er- 
ror indicators, which distinguish between the hierarchical modelling error and the 
finite element error. Application of a posteriori error estimation is important to 
single-layered as well as multi-layered materials, such as metallic thermal protection 
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systems and composite laminates. One of the goals of this research is to investigate 
a posteriori error techniques applied to multi-layered materials. 

5.2 Conclusions 

It is shown that a residual based a posteriori error estimate method can be effectively 
applied to a two-dimensional orthotropic material under steady-state conditions. The 
equilibrated flux method provides more accurate results for the error estimates than 
the average flux method. Although somewhat more computationally expensive than 
the average flux method, the accuracy of the results from using the equilibrated flux 
method are worth the additional effort. The error estimates for this research are 
determined by performing the analysis on a two-dimensional domain, and while the 
procedure for three-dimensional domains may be more complex to implement, com- 
parable accuracy is expected. The directional error indicators are shown to provide 
some insight to separating the error created by the hierarchical model and the finite 
element method, which will be useful in developing an adaptive scheme. Finally, it 
is determined that it is necessary to modify the application of the equilibrated flux 
method to multi-layered materials by using a discontinuous approximate boundary 
flux to obtain accurate results. One such modification is developed and investigated, 
and the error estimates are shown to predict the error in the solution reasonably well. 

5.3 Recommendations 

The research presented in this work may be extended in several different ways. The 
application of this research was limited to two-dimensional steady-state problems. It 
is recommended that future work include three-dimensional work with hierarchical 
modelling in the through-thickness direction. Additionally, extension to transient 
analyses would be the next step. Development of an adaptive scheme using the er- 
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ror estimates discussed in this research should be investigated. The error estimates 
provide the information needed to determine when and where a more accurate ap- 
proximation is needed, and the adaptive scheme can use the error estimates for au- 
tomatic mesh refinement, p-enrichment, or a combination of both. Finally, for the 
multi-layered problems, other methods of approximating the heat flux on the ele- 
ment boundary to provide more accurate error estimates should be investigated. One 
suggested approach is equilibrating the flux on each layer of the element. 
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